Witt et al. EURASIP Journal on Bioinformatics and Systems Biology 201 1 , 201 1 :3 
http://bsb.eurasipjournals.com/content/201 1/1/3 



0 EURASIP Journal on 

Bioinformatics and Systems Biology 

a SpringerOpen Journal 



RESEARCH Open Access 



Modeling time delay in the NFkB signaling 
pathway following low dose IL-1 stimulation 

Johannes Witt 1 *, Sandra Barisic 2 , Oliver Sawodny 1 , Michael Ederer 1 , Dagmar Kulms 2 and Thomas Sauter 1,3 



Abstract 

Stimulation of human epithelial cells with IL-1 (10 ng/ml) + UVB radiation results in sustained NFkB activation 
caused by continuous IKKp phosphorylation. We have recently published a strictly reduced ordinary differential 
equation model elucidating the involved mechanisms. Here, we compare model extensions for low IL-1 doses (0.5 
ng/ml), where delayed IKKp phosphorylation is observed. The extended model including a positive regulatory 
element, most likely auto-ubiquitination of TRAF6, reproduces the observed experimental data most convincingly. 
The extension is shown to be consistent with the original model and contains very sensitive processes which may 
serve as potential intervention targets. 
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Introduction 

The transcription factor NFftB is of central importance in 
inflammation and anti-apoptotic signaling. Upon stimu- 
lation of human epithelial cells with IL-1, NF/tB becomes 
activated due to proteasomal degradation of its cellular 
inhibitor IftBa. This process requires phosphorylation of 
IkBo. by the upstream kinase IKK|3. Since sustained 
NFftB-dependent expression of anti-apoptotic genes con- 
tributes to the maintenance of a range of cancers, its 
activity is tightly regulated and terminated by a negative 
feedback loop, as NFkB promotes I/tBa synthesis. 
Accordingly, various approaches to anti-cancer strategies 
involve inhibition of the NFkB signaling pathway [1]. 

Interestingly, NFkB is converted into a pro-apoptotic 
factor upon stimulation with IL-1 + UVB. The persistence 
of this effect is ensured by sustained NFkB activity [2] 
caused by sustained phosphorylation of IKKP resulting in 
instant phosphorylation and proteasomal degradation of 
newly synthesized IftBa. Chronic IKKP phosphorylation, 
in turn, is due to UVB-induced inhibition of the responsi- 
ble phosphatase PP2Ac [2]. We investigated the details of 
these processes using a systems biological approach, lead- 
ing to the following ordinary differential equation model 
of IKKP phosphorylation and dephosphorylation [3]: 
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[ILR] = -il(t) • fe a [ILR] 
[ILRc] = il(t) • fe a [ILR] - fei[ILRc] 
[IKKp] = fe p [kin][IKK] - fe dp [PP2A][IKKp] 
[PP2A] = -uv(t) • fe uv [PP2A] 



(1) 



with [ILR](0) = 1, [ILRc](0) = 0, [IKKp](0) = 0, [PP2A](0) 
= 1, [IKK] = 1 - [IKKp]. By western blot analysis, we mea- 
sure IKKp obs = IKKp scale lKK. The factor scale_IKK 
describes the unknown ratio between the strength of the 
IKKp band on the gel and the concentration of IKKp. 

The system variables describe the normalized concen- 
trations of IL-1 receptor [ILR], IL-1 receptor complex 
[ILRc], phosphorylated and unphosphorylated IKKP 
([IKKp] and [IKK]), and PP2Ac [PP2A], the inputs il(t) 
and uv(t) describe IL-1 concentration and UVB radiation. 
Due to the normalization, all kinetic parameters of this 
model are given in s' 1 except for k a (nM" 1 s" 1 ), since il(t) 
is given in nM. 

In the original model [3] the effects of the signaling 
cascade are considered to be negligible, so that the 
IKKP kinase [kin] is assumed to have the same time 
course as [ILRc], i.e., [kin] = [ILRc]. Despite this simpli- 
fying assumption the model works well for high IL-1 
doses (10 ng/ml, corresponding to 0.588 nM), indicating 
that the simplification is justified in this dose range. 

While IKKP is very rapidly phosphorylated upon 
10 ng/ml IL-1 stimulation, delayed IKKp phosphoryla- 
tion could be observed upon 0.5 ng/ml or 0.029 nM 
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IL-1 stimulation (Figure 1). This slowly increasing phos- 
phorylation activity is only insufficiently reproduced by 
the original model, so that a more detailed model of the 
signaling cascade appears appropriate for low IL-1 
doses. However, the mechanism causing the signal delay 
is unknown to date. In principle, various reasons for sig- 
nal delay are conceivable. Among the most prominent 
examples are double phosphorylation as occurring, e.g., 
in the MAPK cascade [4], mechanisms with irreversible 
inhibitors and positive feedback mechanisms. Here we 
will investigate these three mechanisms by extending 
the original model by each of them separately to predict 
the most probable explanation for the delay. 

Modeling 

Each of the three potential mechanisms delaying the sig- 
naling cascade following low doses of IL-1 is modeled as 
a separate building block (see Figure 2). It is then used 
to calculate [kin] as a function of [ILRc], where [kin] is 
[Tu] or [Ypp] or [Xa], depending on the potential 
mechanism (see below). 

The positive feedback mechanism describes the auto- 
catalyzed activation of a protein T to the activated form 
Tu. This mechanism is given as 

[f u] = fe u [ILRc] [T] + fe aut [ILRc] [T] [Tu] - fe du [Tu] (2) 

with [Tu](0) = 0 and [T] = 1 - [Tu]. 

In the double phosphorylation mechanism, the active 
form Ypp of a protein Y is reached after two sequen- 
tial modifications (e.g., phosphorylations) catalyzed by 
ILRc: 

[Yp] = feypUILRcllY] -feydpdYp] -feypallLRcHYp] + fe ydp2 [Ypp] 
[Ypp] = k yp2 [ILRc][Yp]-fe yd p2[Ypp] 

with [Y] = 1 - [Yp] - [Ypp], [Yp](0) = 0, [Ypp](0) = 0. 
In the irreversible inhibitor mechanism, an irreversible 
inhibitor I blocks the activated form Xa of some protein 
X within the signaling cascade. The mechanism is 
implemented as 

[X] = -fe act [ILRc][X]+fe dact [Xa] 

[Xa] = fe act [ILRc][X] - fe dact [Xa] - fe inh [Xa][I] ( 4 ) 
[I] = -feinh[Xa][I] 

where [X](0) = 1, [Xa](0) = 0, [I](0) = I„. Note that [I] 
can also be obtained from the mass conservation 
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Figure 1 IKKp phosphorylation is delayed following stimulation 
with 0.5 ng/ml IL-1 (adapted from [3]). 



relation [I] = [X] + [Xa] + I 0 - 1. All concentrations 
except [I] are scaled without loss of generality so that all 
non-zero initial concentrations except I 0 are 1 (see the 
analogous argumentation in Additional File 1 of [3]). 
The exception of [I] is due to the term Ar inh [Xa] [I] 
which occurs in the differential equations for both [Xa] 
and [I] Figure 2 shows a schematic representation of the 
original model and the three extensions. 

While no specific irreversible inhibitor mechanism is 
known within the IL-1 signaling cascade, the double 
phosphorylation mechanism might, e.g., correspond to 
IRAK-1 or IRAK-4 phosphorylation [5], A promising 
candidate for the positive feedback mechanism in the 
IL-1 signaling cascade is TRAF6: The E3-ligase TRAF6 
is a key mediator of IL-1 induced signaling because its 
auto-ubiquitination is required to activate IKKP [6]. 
Hence, the positive feedback mechanism might describe 
TRAF6 auto-ubiquitination. In the following, we there- 
fore refer to this model extension as the TRAF exten- 
sion. Note that, while no other comparable mechanism 
is known to date within the IL-1 signalling cascade, the 
TRAF module might also represent a different, yet 
unknown positively auto-regulated mechanism within 
the signalling cascade. 

Experimental data for 0.5 ng/ml IL-1 stimulation were 
obtained from three western blot experiments measur- 
ing IKI<P Serl77/181 phosphorylation in human epithe- 
lial carcinoma cell line KB following IL-1 stimulation 
with and without UVB (300 J/m ) costimulation until 
120 min post-stimulation as described in [3]. They were 
subsequently extended with data from two western blots 
covering the first 30 min after IL-1 stimulation without 
UVB. Data for the 10 ng/ml IL-1 ± UVB stimulation are 
taken from [3]. Maximal IKKp phosphorylation was 
clearly lower for low (0.5 ng/ml) than for high (10 ng/ 
ml) IL-1 stimulation, so that the maximal value could 
not be used as a reference value for scaling. Therefore, a 
separate scaling factor was used for IKI<P following 0.5 
and 10 ng/ml IL-1 stimulation. 

Modeling and analysis was performed with the Matlab 
(The Math Works) based software tool PottersWheel 
(http://www.potterswheel.de freely available for aca- 
demic use) [7], analogously to the procedures described 
in [3]. Particularly, the % value was used as objective 
function, with 

x 2 (fl)=EE (y *- y( 'r e))1 

i=l j=l l > 

where M is the number of stimulations (i.e., M = 4), 
Ni is the number of data points under stimulation i, y« 
is data point / under stimulation i with standard devia- 
tion Oy and y{tij, 9) is the simulated value at time point /' 
under stimulation i for the parameter vector 9. 
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Figure 2 Graphical representation of the four model variants. Schematic representation of the original IKKp phosphorylation model (a) with 
the extensions positive feedback/TRAF (b), double phosphorylation (c), and irreversible inhibitor (d). Ellipses denote model inputs, boxes with 
solid lines denote state variables, boxes with dotted lines denote species calculated from mass conservation relations. Normal arrows denote 
mass flows, circles at an arrow tail signify that the reactant participates catalytically in the reaction. 



The upper bound of the scaling factor for low IL-1 sti- 
mulations was set to 4, as obtained from experimental 
data (not shown). Upper bounds for the remaining para- 
meter values were derived as described in [3], lower 
bounds were not required. 

Results and Discussion 

Comparison of the different model extensions 

Fitting of the model with its different extensions to the 
experimental data yields better % values than the origi- 
nal model for all extensions (Table 1). Comparison of 
the parameter values with the parameter values of the 
original model showed that the parameters downstream 
of the model extension, i.e., k uv , k P: and kd p , remain basi- 
cally unchanged (Table 2). 

Table 1 % 2 values and AlCc of the model and its 
extensions fitted to the four experimental data sets 



Model 


% 2 value 


AlCc 


Original (a) 


45.8 


124.9 


TRAF (b) 


23.1 


113.7 


Double phosphorylation (c) 


28.3 


123.5 


rreversible inhibitor (d) 


27.0 


122.2 



While all model extensions are able to reproduce the 
delay (Figure 3; Figures SI and S2 in Additional file 1), 
the TRAF extension is the model structure with the 
least degrees of freedom and lowest % value. Conse- 
quently, it also yields the best (lowest) value of the Cor- 
rected Akaike Information Criterion, AlCc, which 
compares the goodness-of-fit of models with a different 
number of parameters (Table 1). In contrast, the other 
extensions only perform about as well as the original 
model in terms of AlCc. 

The TRAF extension convincingly reproduces the 
experimental data and especially the signal delay (Figure 
3), and a biological counterpart to this mechanism exists 
within the IL-1 signaling cascade. In the following, we 
therefore focus on the TRAF extension. All parameter 
values of this extension are given in Table SI in Addi- 
tional file 1. 

Detailed analysis of the TRAF6 model extension 

Fixing of all model parameters also occurring in the ori- 
ginal model to their value reported in [3] yields only 
slightly worse fit quality of % = 24.3. Remarkably, the 
fit quality is still better than the fit quality of the origi- 
nal model [3] {% = 24.5), where only the data sets for 
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Table 2 Values of the kinetic parameters downstream of 
the model extension in the original model [3] and in the 
extended models 



Model 


k p Is" 1 ) 




fcuv (S _1 ) 


Original (a) 


0.095 


7.6e-4 


2.4e-4 


TRAF (b) 


0.095 


9.3e-4 


2.6e-4 


Double phos. (c) 


0.095 


7.5e-4 


2.3e-4 


Irr. inhibitor (d) 


0.095 


7.2e-4 


2.3e-4 



10 ng/ml IL-1 doses ±UVB are considered. As to the 
time course, the most prominent difference of the 
model extension compared to the original model con- 
sists in a longer plateau phase of the maximum concen- 
trations, which is in better agreement to the measured 
data. In summary, the model extension is consistent to 
the model presented in [3] and justifies its additional 
complexity with a considerable enhancement of the fit 



quality and a consistent qualitative behavior for low 
doses. 

Although the parameter fixing performed above sug- 
gests that identifiability in the extended model might be 
worse than in the original one, identifiability analysis 
shows that all parameter values are identifiable from the 
given data: 2000 fits were performed with initial para- 
meter values obtained by randomly perturbing the best 
parameter values by up to four orders of magnitude. 
The standard deviation of the parameter sets of the best 
800 fits is not more than 1% for each parameter. 

Furthermore, we calculated the relative sensitivity 

dx/x 
Sl ~ dpi/Pi 

of x against perturbations in parameter p t following IL- 
1 stimulation for two time-independent characteristics x 




Figure 3 IKKp phosphorylation following different stimulations - experimental data and simulation results of the model with TRAF 
extension. Columns show stimulation with IL-1 or IL-1 + UVB (300 J/m 2 ) radiation, respectively, applied at f = 0 s. Rows indicate whether a high 
(10 ng/ml) or low (0.5 ng/ml) IL-1 dose was applied. A different scaling factor is chosen for the data for 10 and 0.5 ng/ml IL-1 stimulation. 
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of IKKp(f): the IKKp peak amplitude, x = max(IKKp(f)) 
and the mean of IKKp within the longest observed time 
period, 



1 



2.75h 



2.75h 

j IKKp(t) dt 



Oil 



For both sensitivity measures, a high dose dependency 
can be observed (Figure 4): for 10 ng/ml (0.588 nM) IL- 
1, parameter perturbations result in relatively weaker 
changes of the mean IKKp value (|s ; | < 1), and have 
almost no effect on the peak amplitude. In contrast, the 
system is quite sensitive upon 0.5 ng/ml (0.029 nM) IL- 
1 stimulation: for several parameters, perturbations have 
an amplified effect on both peak amplitude and average 
IKKp value (|s ; | > 1). Furthermore, the system is sensi- 
tive over a wide dose range, with sensitivities of up to 
±5 (Figure 4). 



♦ • 



sensitivity 

of peak 
amplitude 

6r 

0.5 ng/ml /*. 

4 l 
2 
0 
-2 
-4 
-6 

b 

sensitivity 
of mean 
cone. 

6r 



V. 



10 ng/ml 

I 



0.02 0.05 /0.10 0.20 

r— 



IL-1 



0.50 1.00 (nM) 



kaut ka ki kdu 



4 

2 
0 
-2 
-4 
-6 



0.5 ng/ml 

1 ^ 
T .♦ 



10 ng/ml 

I 



IL-1 



0.02 



0.05 



0.20 


tr.m-- 


1.00 


kaut 


ka ki 


kdu 





Figure 4 Sensitivity analysis. The relative sensitivity of the (a) 
peak and (b) mean IKKp concentration to perturbations in the 
parameters /r a , /c„ k M , and k du is highly dependent on the IL-1 dose. 
Sensitivities were calculated and plotted using Mathematica 8.0. 



Taken together, the positive feedback regulation pre- 
sumably mediated via TRAF6 auto-ubiquitination [6] 
represents a consistent, well identifiable extension of the 
IKKP phosphorylation model that reproduces the time 
delay at low IL-1 stimulations and is in accordance with 
literature findings. 

On the other hand, it should be noted that determina- 
tion of low levels of phosphorylated IKKP is an experi- 
mentally challenging task, so that the values for 0.5 ng/ 
ml IL-1 stimulation should be regarded as somewhat 
semiquantitative measurements despite the partially low 
standard deviation. Against this background, though the 
postulated model consistently explains the experimental 
data one should neither completely reject the other two 
investigated mechanisms, which also perform well, nor 
lose sight of other potential delay mechanisms. Never- 
theless, our results may provide an indication that the 
reported TRAF6 auto-ubiquitination [6] can indeed be 
interpreted as a positive feedback loop. 

The very good accordance of the relatively small 
model with the experimental data together with the 
reported findings [6] make TRAF6 auto-ubiquitination a 
promising candidate for further research. Several system 
parameters could be characterized as very sensitive at 
low IL-1 doses. The respective processes therefore qua- 
lify as potential intervention targets in cancer therapy. 
The systems biological approach in combination with 
the necessary experimental validation can help to 
further elucidate important molecular mechanisms. 

Additional material 



Additional file 1: Figure SI. IKK(3 phosphorylation following 
different stimulations - experimental data and simulation results of 
the model with buffer extension. Columns show stimulation with IL-1 
or IL-1 + UVB (300 J/m 2 ) radiation, respectively, applied at f = 0 h. Rows 
indicate whether a high (10 ng/ml) or low (0.5 ng/ml) IL-1 dose was 
applied. A different scaling factor is chosen for the data for 10 and 0.5 
ng/ml IL-1 stimulation. Figure S2. IKKp phosphorylation following 
different stimulations - experimental data and simulation results of the 
model with double phosphorylation extension. Columns show 
stimulation with IL-1 or IL-1 + UVB (300 J/m 2 ) radiation, respectively, 
applied at f = 0 h. Rows indicate whether a high (10 ng/ml) or low (0.5 
ng/ml) IL-1 dose was applied. A different scaling factor is chosen for the 
data for 10 and 0.5 ng/ml IL-1 stimulation. Table SI. Parameter values of 
the model with TRAF extension Supplemental table 
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factor 6. 
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